home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Cream of the Crop 26
/
Cream of the Crop 26.iso
/
os2
/
octa209s.zip
/
octave-2.09
/
libcruft
/
minpack
/
r1mpyq.f
< prev
next >
Wrap
Text File
|
1996-07-19
|
3KB
|
93 lines
SUBROUTINE R1MPYQ(M,N,A,LDA,V,W)
INTEGER M,N,LDA
DOUBLE PRECISION A(LDA,N),V(N),W(N)
C **********
C
C SUBROUTINE R1MPYQ
C
C GIVEN AN M BY N MATRIX A, THIS SUBROUTINE COMPUTES A*Q WHERE
C Q IS THE PRODUCT OF 2*(N - 1) TRANSFORMATIONS
C
C GV(N-1)*...*GV(1)*GW(1)*...*GW(N-1)
C
C AND GV(I), GW(I) ARE GIVENS ROTATIONS IN THE (I,N) PLANE WHICH
C ELIMINATE ELEMENTS IN THE I-TH AND N-TH PLANES, RESPECTIVELY.
C Q ITSELF IS NOT GIVEN, RATHER THE INFORMATION TO RECOVER THE
C GV, GW ROTATIONS IS SUPPLIED.
C
C THE SUBROUTINE STATEMENT IS
C
C SUBROUTINE R1MPYQ(M,N,A,LDA,V,W)
C
C WHERE
C
C M IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C OF ROWS OF A.
C
C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C OF COLUMNS OF A.
C
C A IS AN M BY N ARRAY. ON INPUT A MUST CONTAIN THE MATRIX
C TO BE POSTMULTIPLIED BY THE ORTHOGONAL MATRIX Q
C DESCRIBED ABOVE. ON OUTPUT A*Q HAS REPLACED A.
C
C LDA IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M
C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY A.
C
C V IS AN INPUT ARRAY OF LENGTH N. V(I) MUST CONTAIN THE
C INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GV(I)
C DESCRIBED ABOVE.
C
C W IS AN INPUT ARRAY OF LENGTH N. W(I) MUST CONTAIN THE
C INFORMATION NECESSARY TO RECOVER THE GIVENS ROTATION GW(I)
C DESCRIBED ABOVE.
C
C SUBROUTINES CALLED
C
C FORTRAN-SUPPLIED ... DABS,DSQRT
C
C MINPACK. VERSION OF DECEMBER 1978.
C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
C
C **********
INTEGER I,J,NMJ,NM1
DOUBLE PRECISION COS,ONE,SIN,TEMP
DATA ONE /1.0D0/
C
C APPLY THE FIRST SET OF GIVENS ROTATIONS TO A.
C
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 50
DO 20 NMJ = 1, NM1
J = N - NMJ
IF (DABS(V(J)) .GT. ONE) COS = ONE/V(J)
IF (DABS(V(J)) .GT. ONE) SIN = DSQRT(ONE-COS**2)
IF (DABS(V(J)) .LE. ONE) SIN = V(J)
IF (DABS(V(J)) .LE. ONE) COS = DSQRT(ONE-SIN**2)
DO 10 I = 1, M
TEMP = COS*A(I,J) - SIN*A(I,N)
A(I,N) = SIN*A(I,J) + COS*A(I,N)
A(I,J) = TEMP
10 CONTINUE
20 CONTINUE
C
C APPLY THE SECOND SET OF GIVENS ROTATIONS TO A.
C
DO 40 J = 1, NM1
IF (DABS(W(J)) .GT. ONE) COS = ONE/W(J)
IF (DABS(W(J)) .GT. ONE) SIN = DSQRT(ONE-COS**2)
IF (DABS(W(J)) .LE. ONE) SIN = W(J)
IF (DABS(W(J)) .LE. ONE) COS = DSQRT(ONE-SIN**2)
DO 30 I = 1, M
TEMP = COS*A(I,J) + SIN*A(I,N)
A(I,N) = -SIN*A(I,J) + COS*A(I,N)
A(I,J) = TEMP
30 CONTINUE
40 CONTINUE
50 CONTINUE
RETURN
C
C LAST CARD OF SUBROUTINE R1MPYQ.
C
END